function [x, resNE, k, info] = cgls(Aname ,shift, b, m, n, kmax, tol, prnt)
% Solves a symmetric system
%
% :math:`(A^T A + shift I)x = A^T b` or :math:`N x = A^T b`,
%
% where `A` is a general `m` by `n` matrix and shift may be positive or negative.
% The method should be stable if :math:`N = (A^T A + shift I)` is positive definite.
% It MAY be unstable otherwise.
%
% 'gemat' specifies an M-file gemat.m (say), such that:
%
%    * `y = gemat(0, x, m, n)` identifies gemat (if you wish),
%    * `y = gemat(1, x, m, n)` gives :math:`y = A x`,
%    * `y = gemat(2, x, m, n)` gives :math:`y = A^T x`.
%
% The M-file must not be called Aname.m!
%
% USAGE:
%
%    [x, resNE, k, info] = cgls( Aname ,shift, b, m, n, kmax, tol, prnt)
%
% INPUTS:
%    Aname:    name of file
%    shift:    if 0 then `cgls` is Hestenes and Stiefel's specialized form
%    b:        used in the formulas in description
%    m, n:     dimensions of the matrix
%    kmax:     maximum number of iterations.
%    tol:      desired relative residual size, :math:`norm(rNE)/norm(A^T b)`,
%              where :math:`rNE = A^T b - N x`.
%    prnt:     1 gives an iteration log, 0 suppresses it.
%
% OUTPUTS:
%    x:        from the formulas in description
%    resNE:    relative residual for normal equations: :math:`norm(A^T b - Nx)/norm(A^T b)`.
%    k:        final number of iterations performed.
%    info:     gives information on the result:
%
%                * 1 if required accuracy was achieved,
%                * 2 if the limit kmax was reached,
%                * 3 if N seems to be singular or indefinite.
%                * 4 if instability seems likely.  (`N` indefinite and `normx` decreased.)
%
% NOTE:
%
%    When :math:`shift = 0`, `cgls` is Hestenes and Stiefel's specialized form of the
%    conjugate-gradient method for least-squares problems. The general shift
%    is a simple modification.
%
% .. Author: - 01 Sep 1999: First version Per Christian Hansen (DTU) and Michael Saunders (visiting DTU).

x = zeros(n,1);
q = feval(Aname,0,x,m,n);
% Let Aname identify itself

% Initialize
   r    = b;
   s    = feval(Aname,2,b,m,n);                     % s = A'b
   p    = s;
   norms0 = norm(s);     gamma = norms0^2;
   xmax = 0;             normx  = 0;
   k    = 0;             info   = 0;

   if prnt
      head = '    k       x(1)             x(n)           normx        resNE';
      form = '%5.0f %16.10g %16.10g %9.2g %12.5g';
      disp('  ');   disp(head);
      disp( sprintf(form, k,x(1),x(n),normx,1) )
   end

   indefinite = 0;
   unstable   = 0;

%---------------------------------------------------------------------------
%% Main loop
%---------------------------------------------------------------------------
   while (k < kmax) & (info == 0)

      k     = k+1;
      q     = feval(Aname,1,p,m,n);                % q = A p

      delta = norm(q)^2  +  shift*norm(p)^2;
      if delta <= 0, indefinite = 1;   end
      if delta == 0, delta      = eps; end
      alpha = gamma / delta;

      x     = x + alpha*p;
      r     = r - alpha*q;
      s     = feval(Aname,2,r,m,n)  -  shift*x;    % s = A'r - shift x

      norms = norm(s);
      gamma1= gamma;
      gamma = norms^2;
      beta  = gamma / gamma1;
      p     = s + beta*p;

   %% Convergence
      normx = norm(x);
      xmax  = max( xmax, normx );
      info  = (norms <= norms0 * tol) | (normx * tol >= 1);

   %% Output
      resNE = norms / norms0;
      if prnt, disp( sprintf(form, k,x(1),x(n),normx,resNE) ); end
   end %while

shrink = normx/xmax;
if     k == kmax,       info = 2; end
if    indefinite,       info = 3; end
if shrink <= sqrt(tol), info = 4; end
return;
